#front matter
rm(list=ls())
library(foreign)
library(lattice)
library(utils)
library(Hmisc)
library(xtable)
library(mice)
library(pscl)

##set the following to YOUR working directory##
setwd('/Volumes/MONOGAN/psARE/')

#read data
data.0<-read.dta('palsStrictPartic.dta')
#data.0<-read.csv('palsStrictPartic.csv')

#drop aggregate indicators of participation
data<-subset(data.0,select=-c(participation))

#multiple random imputation
imp<-mice(data)

#voting
mods.vote<-glm.mids(vote~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#sign
mods.sign<-glm.mids(po_sign~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#work
mods.work<-glm.mids(po_work~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#contact
mods.contact<-glm.mids(po_contact~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#persuade
mods.persuade<-glm.mids(po_persuade~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#persuade
mods.persuade<-glm.mids(po_persuade~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#attend
mods.attend<-glm.mids(po_attend~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#gave
mods.gave<-glm.mids(po_gave~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#workrel
mods.workrel<-glm.mids(po_workrel~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#demonstrate
mods.demonstrate<-glm.mids(po_demonstrate~cost+hr_age+I((hr_age)^2)+education+dm_income+hr_gender+black+hispanic+other+strength+recruitment+south+attendance+homeowner+wj_unemployed+civic, data=imp, family=binomial(link = "logit"))

#Equation for Social Participation
mods.civic<-lm.mids(civic~cost+hr_age+I((hr_age)^2)+education+hr_gender+black+hispanic+other+strength+recruitment+south+homeowner+wj_unemployed, data=imp)

#Equation for Income
mods.income<-lm.mids(dm_income~cost+hr_age+I((hr_age)^2)+education+hr_gender+black+hispanic+other+strength+recruitment+south+homeowner+wj_unemployed, data=imp)

#Equation for Church Participation
mods.attendance<-lm.mids(attendance~cost+hr_age+I((hr_age)^2)+education+hr_gender+black+hispanic+other+strength+recruitment+south+homeowner+wj_unemployed, data=imp)

#Results
sink('sepLogitResults.txt')
print("vote")
xtable(summary(pool(mods.vote))[,-c(6:11)],digits=4)
print("sign")
xtable(summary(pool(mods.sign))[,-c(6:11)],digits=4)
print("work")
xtable(summary(pool(mods.work))[,-c(6:11)],digits=4)
print("contact")
xtable(summary(pool(mods.contact))[,-c(6:11)],digits=4)
print("persuade")
xtable(summary(pool(mods.persuade))[,-c(6:11)],digits=4)
print("attend")
xtable(summary(pool(mods.attend))[,-c(6:11)],digits=4)
print("gave")
xtable(summary(pool(mods.gave))[,-c(6:11)],digits=4)
print("workrel")
xtable(summary(pool(mods.workrel))[,-c(6:11)],digits=4)
print("demonstrate")
xtable(summary(pool(mods.demonstrate))[,-c(6:11)],digits=4)
print("civic")
xtable(summary(pool(mods.civic))[,-c(6:11)],digits=4)
print("income")
xtable(summary(pool(mods.income))[,-c(6:11)],digits=4)
print("church attendance")
xtable(summary(pool(mods.attendance))[,-c(6:11)],digits=4)
sink()

#Create Figure of Causal Model
png("causal.png")
plot(x=c(0,1),y=c(0,1),type='n',axes=F,xlab="",ylab="")
polygon(x=c(0,0,.2,.2),y=c(.4,.6,.6,.4))
text(x=.1,y=.5,"strictness")
polygon(x=c(1,1,.8,.8),y=c(.4,.6,.6,.4))
text(x=.9,y=.5,"participation")
polygon(x=c(.35,.35,.65,.65),y=c(.8,1,1,.8))
text(x=.5,y=.9,"civic engagement")
polygon(x=c(.4,.4,.6,.6),y=c(0,.2,.2,0))
text(x=.5,y=.1,"income")
polygon(x=c(.35,.35,.65,.65),y=c(.3,.5,.5,.3))
text(x=.5,y=.4,"religious attendance")
arrows(.2,.55,.8,.55,length=.1)
arrows(.2,.55,.35,.9,length=.1)
arrows(.2,.55,.35,.4,length=.1)
arrows(.2,.55,.4,.1,length=.1)
arrows(.65,.9,.8,.6,length=.1)
arrows(.65,.4,.8,.5,length=.1)
arrows(.6,.1,.8,.4,length=.1)
dev.off()

